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Abstract 

In this paper we present a general framework that allows one to study discretiza- 
tion of certain dynamical systems. This generalizes earlier work on discretization 
of Lagrangian and Hamiltonian systems on tangent bundles and cotangent bundles 
respectively. In particular we show how to obtain a large class of discrete algorithms 
using this geometric approach. We give new geometric insight into the Newmark 
model for example and we give a direct discrete formulation of the Hamilton-Jacobi 
method. Moreover we extend these ideas to deriving a discrete version of the max- 
imum principle for smooth optimal control problems. 

We define discrete variational principles that are the discrete counterpart of 
known variational principles. For dynamical systems, we introduce principles of 
critical action on both the tangent bundle and the cotangent bundle. These two 
principles are equivalent and allow one to recover most of the classical symplec- 
tic algorithms. In addition, we prove that by increasing the dimensionality of the 
dynamical system (with time playing the role of a generalized coordinate), we are 
able to add conservation of energy to any (symplectic) algorithms derived within 
this framework. We also identify a class of coordinate transformations that leave 
the variational principles presented in this paper invariant and develop a discrete 
Hamilton-Jacobi theory. This theory allows us to show that the energy error in the 
(symplectic) integration of a dynamical system is invariant under discrete canonical 
transformations. Finally, for optimal control problems we develop a discrete maxi- 
mum principle that yields discrete necessary conditions for optimality. These con- 
ditions are in agreement with the usual conditions obtained from Pontryagin max- 
imum principle. We illustrate our approach with an example of a sub-Riemannian 
optimal control problem as well as simulations that motivate the use of symplec- 
tic integrators to compute the generating functions for the phase flow canonical 
transformation . 
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1 Introduction 

Standard methods (called numerical integrators) for simulating motion take 
an initial condition and move objects in the direction specified by the dif- 
ferential equations. These methods do not directly satisfy the physical con- 
servation laws associated with the system. An alternative approach to inte- 
gration, the theory of geometric integrators[27,7], has been developed over 
the last two decades. These integrators strictly obey some of these physical 
laws, and take their name from the law they preserve. For instance, the class 
of energy-momentum integrators conserves energy and momenta associated 
with ignorable coordinates. Another class of geometric integrators is the class 
of symplectic integrators which preserves the symplectic structure. This last 
class is of particular interest when studying Hamiltonian and Lagrangian sys- 
tems since the symplectic structure plays a crucial role in these systems [3, 1,2]. 
The work done by Wisdom [36, 37] on the n-body problem perfectly illustrates 
the benefits of such integrators. 

At first, symplectic integrators were derived mostly as a subclass of Runge- 
Kutta algorithms for which the Runge-Kutta coefficients satisfy specific rela- 
tionships [31]. Such a methodology, though very systematic, does not provide 
much physical insight and may be limited when we require several laws to 
be conserved. Other methods were developed in the 90's, among which we 
may cite the use of generating functions for the canonical transformation in- 
duced by the phase fiow[8,9] and the use of discrete variational principles. This 
last method "gives a comprehensive and unified view on much of the litera- 
ture on both discrete mechanics as well as integration methods" (Marsden and 
West [26]). Names of variational principles differ in the literature, so we have 
decided to refer to Goldstein[10] in this paper: Hamilton's principle concerns 
Lagrangian systems (i.e., refers to a principle of critical action that involves 
the Lagrangian) whereas the modified Hamilton's principle concerns Hamil- 
tonian systems (i.e., refers to a principle of critical action that involves the 
Hamiltonian). Several versions of the discrete modified Hamilton's principle 
can be found in the literature such as the one developed by Shibberu[32] and 
Wu[38]. For the discrete Hamilton's principle, Moser and Veselov[28] and then 
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Marsden, West and Wcndlandt [26,35] developed a fruitful approach. Also, 
Jalnapurkar, Pekarsky and West [21] developed a variational principle on the 
cotangent bundle based on generating function theory. 

In this paper, we focus on the discrete variational principles introduced by 
Guo, Li and Wu [15,16,17] because the theory they have developed provides 
both a discrete modified Hamilton's principle (DMHP) and a discrete Hamil- 
ton's principle (DHP) that are equivalent. We modify and generalize both 
variational principles they introduce by changing the time discretization so 
that a suitable analogue of the continuous boundary conditions may be en- 
forced. These boundary conditions arc crucial for the analysis of optimal con- 
trol problems and play a fundamental role in dynamics. Our approach not 
only allows us to obtain a large class of discrete algorithms but it also gives 
new geometric insight into the Newmark model [29] . Most importantly, using 
our improved version of the discrete variational principles introduced by Guo 
et al., we develop a discrete Hamilton- Jacobi theory that yields new results 
on symplectic integrators. 

In the first part of this paper (sections 2, 3 and 4), we present a discrete 
Hamilton's principle on the tangent bundle and a discrete modified Hamil- 
ton's principle on the cotangent bundle (section 2), we discuss the differences 
with other works on variational integrators (section 3) and show that we are 
able to recover classical symplectic schemes (section 4). The second part (sec- 
tions 5 and 6) is devoted to issues related to energy conservation and energy 
error. We first show that by considering time as a generalized coordinate we 
can ensure energy conservation (section 5). Then we introduce the framework 
for discrete symplectic geometry and the notion of discrete canonical transfor- 
mations. We obtain a discrete Hamilton-Jacobi theory that allows us to show 
that the energy error in the symplectic integration of a dynamical system is 
invariant under discrete canonical transformations (section 6). Finally, in the 
last part (section 7) we develop a discrete maximum principle that yields dis- 
crete necessary conditions for optimality. These conditions are in agreement 
with the usual conditions obtained from Pontryagin maximum principle and 
define symplectic algorithms that solve the optimal control problem. 

In each part, we illustrate some of the ideas with simulations. In particular 
we show in the first part that symplectic methods allow one to recover the 
generating function from the phase fiow while standard numerical integrators 
fail because they do not enforce the necessary exactness condition. The exam- 
ples presented are the simple harmonic oscillator and a nonintegrable system 
describing a particle orbiting an oblate body. In the second part we look at 
the energy error in the integration of the equations of motion of a particle in 
a double well potential using a set of coordinates and its transform under dis- 
crete symplectic map. In the last part, we use the discrete maximum principle 
to study the Heisenberg optimal control problem. 
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2 Discrete principles of critical action: DMHP and DHP 



In this section, we develop a modified version of both variational principles 
introduced by Guo, Li and Wu [15,16,17] and present the geometry associated 
with them. 



2.1 Discrete geometry 

Consider a discretization of the time t into n instants T = {{tk)ke[i,n]}- Here 
tk+i — tk may not be equal to — t^-i but for sake of simplicity we assume in 
the following that tk+i — tk = t & [l,n]. The configuration space at tk, is 
the n-dimensional manifold Mk and — \J Mk is the configuration space on 
T. Define a discrete time derivative operator on T. Note that may not 
verify the usual Leibnitz law but a modified one. For instance, if we choose 
A^ to be the forward difference operator on TT: 

AUitk) := -ilitk + r) - q{tk)) = A,g, 

T T 

then A^ verifies: 

Ai{f{t)g{t)) = Aifit) ■ g{t) + f{t + r) • Aig{t) . (1) 

2.2 Discrete Hamilton's principle 

Our modified version of the discrete Hamilton's principle derived by Guo, Li 
and Wu [15] is the discrete time counterpart of Hamilton's principle for La- 
grangian systems. Consider a discrete curve of points {qk)kelo,n] and a discrete 
Lagrangian Lfi{qk, Afqf) where A^ is a discrete time derivative operator and 
qk is a function of (g^, ^fe+i)- 

Definition 1 (Discrete Hamilton's principle) Trajectories of the discrete 
Lagrangian system going from {to, qo) to {tn, qn) correspond to critical points 
of the discrete action 

n-l 

= E MQi ^Uk)r , (2) 

in the class of discrete curves {qk)k whose ends are (to, qo) (ind {tn, qn) - In other 
words, if we require that the variations of the discrete action he zero for 
any choice of dqf., and 6qo = 6qn = 0, then we obtain discrete Euler- Lagrange 
equations. 
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Note that if wc do not impose tk+i — tk — tk — tk-i, then the discrete action 
would be defined as: 



n-l 



S^^^L,{qi,AUk){tk^i-h), (3) 

k=0 

but the discrete Hamilton's principle would be stated in the same manner ^ . 

To proceed to the derivation of the equations of motion, we need to specify the 
derivative operator, A^. As we will explain below, its definition depends on 
the scheme we consider. We should also mention that our variational principle 
differs from Guo, Li and Wu's since wc consider that the action has only 
finitely many terms and we impose fixed end points. Such a formulation is 
more in agreement with continuous time variational principles and preserves 
the fundamental role played by boundary conditions. For a discussion on this 
topic, we refer to Lanczos [24] section 15. 



2.3 Discrete modified Hamilton's principle 



As in the continuous case, there exists a discrete variational principle on the 
cotangent bundle that is equivalent to the above discrete Hamilton's principle. 



Definition 2 Let be a discrete Lagrangian on TJA and define the discrete 
Legendre transform (or discrete fiber derivative) ¥L : TAi — > T*M. which 
maps the discrete state space TAi to T*Ai by 

{qtAUt)^{qtpt), (4) 

where 



If the discrete fiber derivative is a local isomorphism, is called regular and 
if it is a global isomorphism we say that is hyperregular. 

If Ld is hyperregular, we define the corresponding discrete Hamiltonian func- 
tion on T*M by 

H,{qi pi) = (Pt ^Ut) - Ld{qt ^Ut) , (6) 



^ In this formulation, the t^s are known, so there are no additional variables. 
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where A^gf is defined implicitly as a function of {qf,pf) through equation (5). 
Let be the discrete action summation: 

n— 1 

= E {{Pt ^Ui) - Ha{qUt)) r , (7) 

fc=0 

where r is to be replaced by t^+i — tfc if tk+i —tk 7^ tfc —tk-i- Then the discrete 
principle of least action may be stated as follows: 

Definition 3 (Discrete modified Hamilton's principle) Trajectories of the 
discrete Hamiltonian system going from (to, 5o) to {tn, Qn) correspond to 
critical points of the discrete action in the class of discrete curves {q'^,Pk) 
whose ends are {to,qo) and {tn,qn)- 

Again, for deriving the equations of motion we need to specify the discrete 
derivative operator, and its associated Leibnitz law. It will generally de- 
pend upon the scheme we consider as we will see through examples later. 



3 Compcirison with other classical veiriational principle 

At this point it is of interest to compare discrete variational principles intro- 
duced in this paper and other classical discrete variational principles. As we 
mentioned above, the discrete variational principles we develop are inspired 
by the work of Guo, Li and Wu [15] and we explained above the key difference 
between our work and this earlier work. We now point out the main differ- 
ences of the work discussed here with that of Marsden and West, based on 
the variational principle introduced by Moser and Veselov. In the following, 
DVPI refers to the discrete variational principle developed by Moser, Veselov, 
Marsden, Wendlandt et al. whereas DVPII denotes the discrete variational 
principles developed by Guo and this paper. 

The first main difference lies in the geometry of both variational principles. 
Whereas the discrete Lagrangian is a functional on Q x Q where Q is the 
configuration space in DVPI, it is a functional on TQ in DVPII. As a con- 
sequence, DVPII has a form more like that of the continuous case but has a 
major drawback: we have to specify the derivative operator and the Leibnitz 
law it verifies in order to derive discrete Euler-Lagrange equation. Such a law 
allows us to perform the discrete counterpart of the integration by parts and 
depends on the scheme we consider. On the other hand, the Euler-Lagrange 
equation obtained by DVPI is scheme independent. One benefit is that these 
equations ensure satisfaction of physical laws such as Noether's theorem for 
any numerical scheme which can be derived from them. 

The next important difference between the two discrete variational principles 
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lies in the role of the Legendre transformation in defining a discrete Hamilto- 
nian function from the discrete Lagrangian. In DVPl, one defines a discrete 
Legendre transform to compute the momenta from the discrete Lagrangian 
function, so one may study the discrete dynamics on both Q x Q and T*Q. 
However, it does not seem possible to define a discrete Hamiltonian function 
from the discrete Lagrangian and develop a DMHP. Given a Hamiltonian sys- 
tem, to derive discrete equations of motion using DVPl one needs to first find 
a continuous Lagrangian function by performing a Legendre transform on the 
continuous Hamiltonian function, then apply DVPl and finally use the discrete 
Legendre transform to study the dynamics on T*Q (see for instance [26] page 
408). While this point may not be of importance when dealing with dynami- 
cal systems, it is crucial if one wants to discretize an optimal control problem, 
where the continuous Hamiltonian function does not have any physical mean- 
ing and the Legendre transformation may not be well-defined (See section 7) . 
DVPH naturally defines a discrete Legendre transform and a DMHP. 

As mentioned in the introduction, people have already introduced DMHPs 
on the cotangent bundle, but, as far as we know, no one has developed an 
approach that allows one to equivalently consider both the Hamiltonian and 
Lagrangian approaches in discrete settings (i.e., a DMHP and a DHP that are 
equivalent for non-degenerate Lagrangian systems). In addition, the DMHPs 
that can be found in the literature do not allow one to recover most of the 
classical schemes. For instance, Shibberu's DMHP focuses on the midpoint 
scheme and Wu developed a different DMHP for each scheme. 

Let us now look at some classical schemes and see how they can be derived 
from DVPH. 



4 Examples 

4-1 Stormer's rule and Newmark methods 

Stormer's scheme is a symplectic algorithm that was first derived for molecular 
dynamics problems. It can be viewed as a Runge-Kutta-Nystrom method in- 
duced by the leap-frog partitioned Runge-Kutta method[31]. The derivation of 
Stormer rule as a variational integrator came later and can be found in [38,35]. 
Guo, Li and Wu [17] recovered this algorithm using their discrete variational 
principles. In the next subsection, we briefiy go through the derivation and 
add to their work the velocity Verlet [34] and Newmark methods [26]. In par- 
ticular, we will show how the conservation of the Lagrangian and symplectic 
two- form is built into DVPH. 
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4.1.1 From the Lagrangian point of view 



We first let = qk and define the discrete Lagrangian by Ld^q^i ^t^k) ~ 
L{qk,Afqk) and the discrete derivative operator as the forward difference = 
At-. A^ satisfies the modified Leibnitz law (1). Discrete equations of motion 
are obtained from discrete Hamilton's principle (definition (1)): 



n-1 



5S^^rJ2SHQk,^rqk) (8) 

k=0 
n-1 

^'Tj2(^i^d.{qk:^Tqk):Sqk) + {D2Ld{qk,Kqk),5Kqk) (9) 

k=Q 
n-1 

^^^{DiLa{qk, Kqk) - A^L>2-t'd(gfe-i, Kqk-i), ^qk) 

k=l 

+Ar{D2Ld{qk-i, Arqk-i), Sqk) 

+T{DiLd{qo, A^qo)Sqo) + TL>2Ld(go, A^qo)5A^qo (10) 

rt— 1 

^T^{DiLd{qk,Arqk) - ArD2Ld{qk-i, Arqk-i),Sqk) - 
k=l 

-{D2Ld{qo, Arqo), Sqo) + T{DiLd{qo, Arqo)Sqo) 

+ {D2Ld{qn-U Ar^n-l), (^?n) , (H) 

where the commutativity of 6 and A,- and the modified Leibnitz law defined 
by equation (1) have been used. 

Discrete Euler-Lagrange equations follow by requiring the variations of the 
action to be zero for any choice of 5qk, k & [l,n — 1] and Sqo — 5qn — 0: 

DiLdiqk, Arqk) - ArD2Ld{qk-i, A^qk-i) = . (12) 

Suppose L{q, q) — ^qMq — V{q), then equation (12) yields Stormer's rule: 

qk+i = 2qk - qk-i + h?M-\-VV{qk)) . (13) 

Consider the one- form ^ 

aL dLd(qk-i, Arqk-i) , i 

— aK:^, — 

and define the Lagrangian two- form uj^ on Tq^ Ai: 



^ Einstein's summation convention is assumed 
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= d^k 

^ a i aA J '^^'i ^ + aA i aA i ^^^^^^ ^ -(^^^ 

dql_idArqi_i dKo^k-id^rqi-i 

Lemma 4 T/ie algorithm defined by Stormer's rule preserves the Lagrangian 
two- form, uj^. 



PROOF. Consider a discrete trajectory {qk)k that verifies equation (13). 
Then we have: 



^qL (dLd{qk,Arqk) dLd{qk-i,A^qk-i)\ ^ 
dSd-r}_^l — A, — — dg. 



fc=i 



dqi dAdqi_, 



9AA "V- 

Since the q^s verify equation (13), and d"^ = 0, equation (15) yields: 

d{ArOk) = , that is, ooj^^^ = ooj: . (16) 
We conclude that is preserved along the discrete trajectory 



As we mentioned earlier, because DVPII acts on the tangent bundle it provides 
results very similar to the continuous case as attested by the form of the 
Lagrangian 2-form. This is to be compared with the Lagrangian two-form 
arising in the continuous case: 



Note that conservation of the Lagrangian two-form is a consequence of us- 
ing the Leibnitz law, and therefore does not depend on the definition of the 
discrete Lagrangian. In the remainder of this section we use different discrete 
Lagrangian functions, but the same Leibnitz law. Thus lemma 4 still applies. 

More generally, we can derive Stormer's rule using 

Ld{qk, Arqk) = \L{qk, A^g^) + (1 - X)L{qk + rA^^fc, A^g^) , 

for any A in M. A particular case of interest is A = | which yields a symmetric 
version of Stormer's rule also called the velocity Verlet method [34]. For this 
value of A, we define the associated discrete momenta using the Legendre 
transform (equation (5)): 
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Pk+i^pi (18) 

^D2Ld{qk,Arqk) (19) 
= MArqk - ^rWV{qk + A^q^) , (20) 

that is: 

qk+i = Qk + rM-Hpk+i + -TVV{qk+i)) . (21) 
Moreover, from equation (12) we obtain: 

^ -VVjqk) - VV{qk+i) ...^ 
Pk+i =Pk + r . (22) 

Equations (21) and (22) define tfie velocity Verlet algorithm. 

We now focus on the Newmark algorithm which is usually written for the 
system L = ^q^Mq - V{q) as a map given by {qk, %) ^ {qk+iAk+i) satisfying 
the implicit relations: 



r2 

qk+i = qk + rqk + y [(1 - 2p)ak + 2(3ak+i\ , (23) 

qk+i = qk + r[{l--f)ak + 'yak+i], (24) 
ak^M-\-VV{qk)), (25) 

where the parameters 7 e [0,1] and P e [0,1]. For 7 = | and any /3 the 
Newmark algorithm can be generated from DVPII as a particular case of the 
Stormer rule where q^ and are chosen as follows: 

(it^Qk- Pr^dk , 

and 

L,{qt AUt) = \q{Mqi V{qt) , 

with V, the modified potential, satisfying VV(g^) = W{qk). Since the deriva- 
tive operator is the same as above, the discrete Hamilton's principle yields 
Stormer 's equation where q^ is replaced by , that is: 

4+1 = ^4 - 4-1 + T'M-\-VV{qt)) . (26) 
Equation (26) simplifies to 

qk+i - 2qk + qk-i = r^(/5afc+2 + (1 - 2P)ak+i + Pak-i) . (27) 

This last equation corresponds to the Newmark algorithm for the case 7 = |. 
Lemma 4 guarantees that the Lagrangian two-form 

u;^^d{D2L,{qi,Ayk)dqU,) 
is preserved along the discrete trajectory. 
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4.1.2 From the Hamiltonian point of view 



The Stormer, velocity Verlet, and Newmark algorithms can also be derived 
using a phase space approach, i.e., the DMHP. For Stormer's rule, the Legendre 
transform yields: 

Pfe+i = MArQk . (28) 
The discrete Hamiltonian function is defined from equation (6): 

Hd{qk,Pk+i) = ^pI+iM- Vfe+i + V{qk) , (29) 

and discrete equations of motion are obtained from the DMHP^ (theorem 
(3)). We skip a few steps in the evaluation of the variations of to finally 
find: 



SS^-S fr|:(pfc+i, A,gfe) - Ha{qk,Pk+i)) (30) 

\ k=0 J 
n-1 

^rj^i^rqk - D2Hd{qk,Pk+i), Spk+i) - {^rPk + DiHd{qk,Pk+i), Sqk) 

k=0 

+ {Pn, Sqn) - {po, Sqo) . (31) 

If we impose the variations of the action to be zero for any {5qk, Spk+i) 
and Sqo — Sqn — 0, we obtain: 



Arqk^Pk+i, (32) 
A,pfc = -Vy(gfe). (33) 

Elimination of the p^s yields Stormer's rule. 

To recover the velocity Vcrlet scheme from the Hamiltonian point of view, one 
needs to solve for ArQk as a function of {qk,Pk+i) in equation (20). Suppose 
this has been done and that A^-g^ = f{qk,Pk+i), then 

Hd(qk,Pk+i) = {Pk+i, f{qk,Pk+i)) - Ld(qk, f(qk,Pk+i)) ■ (34) 

Taking the variation of the action Sd yields the following discrete Hamilton's 
equations: 



Arqk = D2Hd{qk,Pk+i) , (35) 
ArPk = -DiHd{qk,Pk+i) ■ (36) 



^ (lk = 9k and pf = pk+i 
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On the other hand, equation (34) provides the following relationships: 

DiHd{qk,Pk+i) = Dif{qk,Pk+i)iPk+i - D2Ld{qk, f{qk,Pk+i))) 

- DiLd{qkJ{qk,Pk+i)) , (37) 



D2Hd{qk,Pk+i) = ^Uk + D2f{qk,Pk+i){Pk+i - D2Ld{qk, f{qk,Pk+i))) ■ (38) 

Combining equations (35) and (36) together with equations (37) and (38) 
yields the Velocity Verlet algorithm (equations (21) and (22)). 

We now prove that the scheme we obtained is symplectic. As in the Lagrangian 
case, the proof differs from the usual one that consists in computing dpk+i A 
dqk+i, in that it relies on fundamental properties of DVPII and on the use of 
the Leibnitz law. 

Lemma 5 The algorithm defined by equations (35)- (36) is symplectic. 
PROOF. We have: 

dSd ^dirY^ipk+u^rqk) - Hd{qk,Pk+i)] , (39) 

\ k=0 J 

n—1 

^-T^i^rqk - D2Hd{qk:Pk+l):dpk+l) - {A^Pk + DiHd{qk: Pk+l) , dqk) 
k=0 

+A^{pk,dqk). (40) 
Hence, since {qk,Pk) verifies equations (35)-(36) and d^ — 0, we obtain: 

A^{dpkAdqk) = 0. (41) 
The symplectic two-form dpk A dqk is preserved along the trajectory. 



Finally, we can also derive the Newmark methods from the Hamiltonian point 
of view. The Legendre transform yields: 

Pt-^^^^^^^^-MAUi. (42) 

The Newmark algorithm is again a particular case of the Stormer rule where 
(qk,Pk+i) is replaced by {qtpt)- 



^Ut=Pt, (43) 
Ay, = -VV{qt) . (44) 
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Defining from pk as 



allows us to recover the Newmark scheme for 7 = | (equations (23) and (24)). 
Prom the above lemma, we obtain that the symplectic two- form dpf A dq^ is 
preserved along the trajectory. 



4-2 Midpoint rule 



The midpoint rule has been extensively studied and a complete study of its 
properties can be found in the literature. It is a particular case of the Runge- 
Kutta algorithm, but can also be derived as a variational integrator (sec for 
instance [38,32,26]). The derivation of this scheme has been done by Guo, 
Li and Wu [17] for the Hamiltonian point of view. In the next section we 
present the Lagrangian point of view and then recall the Guo, Li and Wu 
main results, the goal of this section being to illustrate the use of DVPII with 
other discretization and discrete derivative operator. 



4.2.1 From the Lagrangian point of view 

Given a Lagrangian L{q, q), define the discrete Lagrangian by: 

L,{qt, AUt) = L{qi AUt) , (45) 

where = ££+i±£fe^ = Rr/2 — R-t/2 where the operator is the 

translation by r. One can readily verify that A^g^ = A^-^fc and that A^ verifies 
the usual Leibnitz law: 

A' Utgt) = Kfk ■ 9i + ft ■ ^gi , (46) 

where /. = /(t.) and gu = git,) are functions of time and = 
Applying the discrete Hamilton's principle yields: 



n-1 



5S^^Tj:SL,{qtAUt) (47) 

fe=0 

n-1 

^Tj:{D,L,{q',, AUil 5qi) + {D,La{qt AUt),5AUt) • (48) 



fe=0 



Prom the Legendre transform (equation (5)), we define the associated momen- 
tum: 

^-^^^^pt^D,LM,^M- (49) 
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Then, equation (48) becomes: 



SS^,=TY.{D,La{ql/\Ut).^4) + ipi^^Ut) (50) 

fc=0 

n— 1 

= TY.{D,L,{qt Ay,) - A^pt, 6qi) + - (po, 5go) • (51) 

fc=0 

If we require the variations of the action to be zero for any choice of k G 
[1, n — 1], and 5qQ = Sqn = 0, we obtain discrete Euler-Lagrange equations for 
the midpoint scheme: 



^D,L,{qlAUi) 

n T ^Qk+i + qk qk+i-qk^ frr,\ 

= DiLd{ , ) , (52) 



Pk+i +Pk d 
-^^-Pk 



^D^Ld{qlAUi) 

n T /^k+i + qk qk+i — qk-. /rQ\ 
= D2Ld{ — - — , — - — ) . (53) 

Lemma 6 The midpoint scheme defines a symplectic algorithm. 



PROOF. The proof proceeds as for the Stormer rule: 



n-l 

dS^ ^tY: {D,Ld{qt Aiqi), dqt) + {pi dAfqt) (54) 

k=0 

= T j:{D,Ld{qt, AUt) - A^pi, dqt) + A%1 dqt) . (55) 

k=0 

Since d"^ — and {qk,Pk) verifies equations (52)-(53), we obtain: 

A'^idpi A dqt) = . (56) 

A straight forward computation shows that Af{dpf A dq"^) = Ar^dp^ A rfg^), 
i.e., the symplectic two-form uj^ — dpk A dqk is preserved along the trajectory. 



4-2.2 From the Hamiltonian point of view 

Let Hd{q'l,pf.) = H{qf,pf.) or equivalently define Hd from Ld via equation (6) 
and let {q^pf) = {3^+^^ P£±i±Pfc) xhen the DMHP (3) yields: 
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h 




(57) 



Pfc+l ~ Pk 
h 




)■ 



(58) 



Lemma 7 The midpoint scheme defines a symplectic algorithm. 



PROOF. The proof is straightforward. We compute d'^S^^ assuming {qk,Pk) 
verifies the above equations of motion. 



To conclude, we have illustrated the use of the discrete variational principles 
(definitions (1) and (3)) and derived discrete equations of motion. One can 
readily verify that both variational principles yield the same discrete equa- 
tions, as in the continuous case. Other schemes can be recovered in the same 
way, and we do not know yet if all classical symplectic algorithms can be 
derived from DVPII. For instance, we have been able to recover the condi- 
tions for the partitioned Rungc-Kutta algorithm to be symplectic from the 
Lagrangian point of view but so far it is not clear to us how it can be done 
using the Hamiltonian approach (definition (3)). 



4-3 Numerical example 

Symplectic integrators are usually used as numerical integrators that preserve 
the qualitative behavior of dynamical systems and are especially valuable for 
long time simulations. However, these are not the only uses of symplectic in- 
tegrators. In this section we present an aspect of symplectic integrators that 
we have not seen pointed out in the literature: we show that they allow one 
to recover the generating functions for the phase fiow canonical transforma- 
tion, whereas numerical integrators do not, even over a short period of time 
(applications of this result can be found in [14]). 

Let us first recall two results from the Hamilton- Jacobi theory. 

Proposition 8 The transformation induced by the phase flow is canonical. 
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Proposition 9 Let (Pi,c<Ji) and (^2,^2) be symplectic manifolds, tt^ : Pi x 
P2 Pi the projection onto P^, i — 1,2, and 

Q = TT^Wi — n2UJ2 ■ (59) 

Then: 

(1) D, is a symplectic form on Pi x P2; 

(2) a map / : Pi — > P2 is symplectic if and only if 1*^0. — 0, where : F/ — > 
Pi X P2 is inclusion and is the graph of f. 

Hence, by the Poincare lemma, if / is canonical there exists a function S such 
that 

i}e = dS , (60) 

where Q = —dO. S is called a generating function. If {q^,Pi) are coordinates on 
Pi and (Q*, Pj) are coordinates on P2, then Ff can be endowed with a chart in 
several ways. For instance, S may appear as a function of {q\ Q^) or of {q\ Pi), 
and so forth depending of the choice of 0. Let 61 = pidq^ and 62 = PidQ\ 
then i*fQ = ijirlOi - 1)7^192 = (tti oifYp^dq'' - [712 o if)*PidQ\ In this case, S 
is a function of {q^, ■ ■ ■ ,q",,Q^, - • ■ , Q"). From 

we conclude, using equation (60) that: 

Suppose that / is the phase flow $, then equation (61) defines a relationship 
between flow and the gradient of the generating function. In particular, if the 
generating function S{q, qo, t) exists and the flow is defined as follows: 

$ : iqo,Po,t) ^ {q{t),p{t)) = {Mqo,Po),MQo,Po)) , (62) 

then, from the local inverse function theorem ^ , there exist two functions 
and S2 such that: 



Po^Si{q,qo,t) , (63) 
p^^2{qo,Si(q,qo,t)) = S2(q,qo,t) . (64) 

From equation (61), we conclude that 5*1 and 5*2 are the gradient of S and 



4 I 



^ I 7^ since we assume that S exists 
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therefore should verify ^ : 



d'^S ^ dSi ^ _ ^/ q t)^ .ggx 
dqodq dq ' ' dqo ' ' dqdqo 

Since symplectic integrators preserve the symplectic two-form, the exactness 
condition (equation (65)) is satisfied whereas it is not using numerical inte- 
grators. 



4.3.1 Harmonic Oscillator 

We start with a trivial example, the harmonic oscillator, because its study 
allows us to introduce techniques and discuss issues that arise in the next more 
sophisticated example. The Hamiltonian function for the harmonic oscillator 
is quadratic: 

HM = ^P' + ^l'- (66) 
It is a linear system so the phase flow is also linear: 



^i{qo,Po) = aii(t)qo + ai2(t)po (67) 
•^2(^0,^0) = 021 (t) go + a22(t)po ■ (68) 

Substituting these expressions into Hamilton's equations and balancing terms 
of the same order yield: 



aii{t) 


= a2i{t)/m 


di2{t) 


= a22(t)/m 




= kaii{t) 


_ a22{t) 


= kai2{t) 



(69) 



In figure 1, wc plot A = ^(g, go, ^) ^^(.1- ^) over the time interval [0, 100] 
using the midpoint scheme with fixed time step, a symplectic Gauss implicit 
Runge-Kutta algorithm of order 8 with fixed time step and a non symplectic 
Runge-Kutta integrator of order 8 to integrate equations (69). We remark 
that only symplectic integrators allow us to recover the generating functions 
because the exactness condition is exactly verified. We point out that even over 
a short time span, numerical integrators fail to satisfy the exactness condition. 



^ Since their exists an open set on which the generating functions are smooth, 



Schwartz's theorem yields 



dqodq dqdqo ' 
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(a) Midpoint scheme with fixed 
time step r = 0.01 



(b) Implicit Gauss Runge-Kutta 
algorithm of order 8 



0.0004 
0.0002 




(c) Explicit Runge-Kutta algo- 
rithm of order 8 

Fig. 1. Exactness condition using 3 different integrators 

4.3.2 Earth orbit 

This example was first encountered by V.M. Guibout and D.J. Scheeres [12,14] 
while studying spacecraft formation flight. Consider an orbital problem about 
the Earth modelled by a non-spherical body (we take into account J2 and J3 
gravity coefficients) . The Hamiltonian of the system is given by 



' ^ (3 I 

y/x'^ + y'^ + z'^ \ 2rQ{x^ + + z^) + + z'^ 

z^ \ 
~2rg(x2 + y2 + ^2)2 (5^2 + ^2 + ^2 - -^3 J , 

where 

GM = 398, 600.4405 km^s'^ , i? = 6, 378.137 km , 
J2 = 1.082626675 • 10"^ , J3 = 2.532436 • lO'^ , 
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and all the variables are normalized (ro is the initial radius of a trajectory): 



GM 



X — > XTq . 

(71) 

Pz 






We choose the nominal trajectory to be a highly eccentric orbit. The initial 
conditions for the nominal trajectory are in normalized units {tq — 7000 km): 

X — 1 , y — 0, z — 0, 

(72) 

= , py^ ^ cos(7r/3) , p^ = sin(7r/3) . 
At the initial time, e = 0.3, i — \ Tad, uj — Q and Q = 0. 

This system is non-integrable and has non-trivial dynamics. The phase flow is 
not known globally but techniques have been developed to evaluate it locally 
(see Guibout and Scheeres [13]). 

Consider a given trajectory called the nominal trajectory (5'°(t),p°(t)), then 
the dynamics of the relative motion of a particle about this trajectory is Hamil- 
tonian and is described by the following Hamiltonian function H^{X^,t): 

E2ri 

(73) 

where is the relative state vector (Agr, Ap). 

In the same way, we expand in Taylor series the phase flow for the relative 
motion, and substitute its expression into Hamilton's equations for H^. When 
studying spacecraft formation flight we often assume that the spacecraft stay 
close to each other and therefore, one may approximate the dynamics of the 
formation by truncating the above Taylor series. Suppose we keep only terms 
of order less that N . Then balancing terms of the same order in Hamilton's 
equations yields a set of ordinary differential equations (the procedure is the 
same as in the harmonic oscillator example but here there are non linear terms 
up to order N). We use the midpoint scheme with fixed-time step (r = 0.01), 
a symplectic Gauss implicit Rungc-Kutta algorithm of order 4 with fixed time 
step (r = 0.01) and MathemafAca® built-in numerical integrator NDSolve^ 
to integrate the flow up to order = 4. Once the Taylor series of the flow is 
known, we find Si and S2 by a series inversion. Then we check the exactness 
conditions defined by equation (65) (there are several terms involved since 
we are dealing with a nonlinear system of dimension 6 ). We find that after 



^ NDSolve switches between a non-stiff Adams method and a stiff Gear method. 
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IOtt units of time, ||^ — ^|| < V-i where rj = 10 using the midpoint 

scheme, 1] = 10~^^ with the symplectic imphcit Gauss Runge-Kutta algorithm 
and 1] = 10^'^ with the buih-in function NDSolve. Again, only symplectic 
algorithms allow us to recover the generating functions. 



5 Energy conservation 

Symplectic integrators do not conserve energy and in general induce bounded 
energy error. There are several works on analyzing the energy error, we refer 
to Hairer and Lubich [19] and Haircr, Lubich and Wanner [20] and references 
therein for more details. In this section, we enhance DVPII so that energy con- 
servation is imposed. By considering the time as a coordinate and by adding 
an independent parameter r, DVPII yields symplectic energy conserving al- 
gorithms. For certain problems, such algorithms may provide better perfor- 
mance^, but the contrary may also happen [18,33]. The method we develop 
in this section is variational and allows us to recover Shibberu's algorithm [32] 
for Hamiltonian systems and is equivalent to the Kane, Marsden and Ortiz 
[23] method for Lagrangian systems. 

5. 1 Generalized variational principles 

5.1.1 Generalized Hamilton's principle 

Let us first recall Hamilton's principle for dynamical systems for which time 
is considered as a generalized coordinate. Such a formulation is typically used 
in relativity where the time coordinate is equivalent to the space coordinates. 

Consider a Lagrangian L(g, q) and define the parametric Lagrangian 

L{q,t,ci,t')^t'L{q,^,t), 

where ' = ^ and r is an independent parameter that parameterizes the tra- 
jectory and the time. Then the generalized Hamilton's principle reads: 

Definition 10 Critical points of jl^ L{q, t)dT in the class of curves (^(t), i(T)) 
with endpoints {qo,to) and {qf,tf) correspond to trajectories of the Lagrangian 
systems going from {qo,to) to {qf,tf). 

To quantify the performance of an algorithm, not only we look at its accuracy 
but we also evaluate its ability to predict the qualitative behavior of the system. 
In that sense, symplectic-energy conserving algorithms may not predict qualitative 
behavior better that symplectic algorithms. 
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The generalized Hamilton's principle yields the following set of equations: 



Replacing the parametric Lagrangian by the Lagrangian of the system simpli- 
fies the above equations to: 



^,dL d dL ^ 
dq dr dq 

These n + 1 equations should be compared to the n equations obtained when 
the trajectory is parameterized by the time: 

Since jp — t'-^, we conclude that the space components of the generalized 
Euler-Lagrange equations (equation (77)) are a multiple by t' of the orig- 
inal Euler-Lagrange equations (equation (78)). Also, their time component 
(equation (76)) is a linear combination of the components of equation (78) 
(the sum of each component multiplied by q') . All n + 1 generalized Euler- 
Lagrange equations are thus consistent with the original equations but there 
is no unique solution because they are satisfied by any parameterization. To 
get a unique solution, it is necessary to add to the generalized Hamilton's 
principle an additional condition fixing the parameterization. As we will see 
in the next section, in discrete settings we do not have this freedom anymore. 
The discrete counter-part of equation (76) corresponds to an energy constraint 
that fully specifies the time parameterization, i.e., the time step. 



5.1.2 Generalized discrete Hamilton's principle (GDHM) 

In contrast with the variational principles introduced in the first part of this 

paper, we do not set the time step, i.e., we let the time act as a variable by 
adding an independent parameter such that tk = t{rk) and r^+i — Tk = r, 
T being a constant, tk is now a coordinate that plays the same role as qk, 
is the extended configuration space {qk,tk), M. — {jM^ and T = {{Tk)k£[i,n]}- 
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Define the modified discrete Lagrangian L^: 

LM, 4, ^Ut Kti) = ^%LMt |||, 4) , (79) 

wliere is the discrete Lagrangian previously defined. In addition, since we 
are interested in conservation of energy we only consider system that are time 
independent. As a consequence, does not depend on time and = 0. 

Definition 11 (Generalized Discrete Hamilton's Principle (GDHP)) 

Critical points of the discrete action 

11—1 

S^^j:L,{qitiAUt,AX)T, (80) 

fc=0 

in the class of discrete curves (g^, t^)^ with endpoints (tq, to, go) (ind (r„, t„, g„) 
correspond to trajectories of the discrete Hamiltonian system going from (to, go) 
to itn,qn): 

Again, to proceed to the derivation of the equations of motion we need to 
specify the derivative operator. 



5.1.3 Generalized discrete modified Hamilton's principle 

Definition 12 Let he a discrete Lagrangian on TM. and define the discrete 
Legendre transform (or discrete fiber derivative) ¥L : TM. — > T*A4 which 
maps the discrete eoctended phase space TAi to T*Ai by 

{qtttAUk:AX)-^{qitipiet), (81) 

where 

, _ dUqitiAjqlAitt) , _ dL.iqltiAjqlAX) 

The Legendre transform as defined by equations (82) is equivalent to the pre- 
vious definition (equation (5)^. Indeed, 

d ^UU 



^' -D,L,{qt,AUt), 



dAUi dAUt 
where — represent the discrete derivative with respect to time. 



If the discrete fiber derivative is a local isomorphism, L^ is called regular and if 
it is a global isomorphism we say that L^ is hyperregular. If L^ is hyperregular. 
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we define the corresponding discrete Hamiltonian function on T*M by 

H,{qttt,pt,et) = {pi Ay,) - UqtAUt) , (83) 

where Afq, is defined implicitly as a function of {qi,pf) through equation (82). 
Hd is related to the previously defined Hamiltonian function by the following 
relationship: 

Haiqlvi) ^ A%Ha{qUi) . (84) 
In addition, we have: ef = —Hd{q^,pf.)j that is, the momentum associated 
with the time is the opposite of the Hamiltonian. 

Let be the discrete action summation: 



n-1 



^rY.{pt,Aiqi) - H,{qipi) 



k=0 
n-1 



■T^{pt,AUi) + {etAiti) 

k=0 



(85) 
(86) 



Before stating the generalized discrete modified Hamilton's principle, we need 
to remark that all the coordinates are not independent since the holonomic 
constraint ef. = —H{q^,pf) holds. There are two ways to handle this situa- 
tion [3], one can either replace ef. by —H{qf,pf) in the action and then take 
the variations or one can use Lagrange multiplier to append the constraint 
ef. + H{q'^,pf.) = to the integral. Therefore we can give two equivalent for- 
mulations of the GDMHP. 

Definition 13 (Generalized discrete modified Hamilton's principle) 

Critical points of the discrete action 

n-1 

Sf-rj:{pt,Aiqt) + {et,Aitt) 

k=0 

in the class of discrete curves {q^,tf,pief) with endpoints {TQ,to,qo) and 
{jri,tn,qn) subject to the constraint ef -|- Hd^qt^Pk) = correspond to discrete 
trajectories of the discrete Hamiltonian system going from {tQ^qo) to {tn,qn)- 

Definition 14 (Generalized discrete modified Hamilton's principle) 

Critical points of the discrete action 

n— 1 

Si = rj^ipt, AUt) - Ha{qt,pt)AX 

k=0 

in the class of discrete curves {qf, tf.,pf) with endpoints (tq, to, go) and (r„, tn, qn) 
correspond to trajectories of the discrete Hamiltonian system going from {to, qo) 

to {tn,qn)- 
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To derive the equations of motion wc need to specify the discrete derivative 
operator, and its associated Leibnitz law. 



5.2 Examples 



5.2.1 Stormer type of algorithm 

5.2.1.1 Lagrangian approach Consider a Lagrangian function L(q,q) 
and define the discrete Lagrangian map trivially by Ld{qk, ^rQk) = L{qk, ^rQk)- 
Discrete equations of motion are obtained from the generalized discrete Hamil- 
ton's principle: 



n-1 



^^d = S ^Ld{qk, tk, A^qk, ^rtk) 



k=0 
n-1 



= T ^ S{ArtkLd{qk, ^^)) 
k=Q ^rT^k 



n-1 



k=0 



where = Ld{qk, x^)- Using the Leibnitz law (equation (1)) and the fixed 
end points constraint, we obtain: 

n— 1 

^tJ2 -^rCk-iStk + (-A,L>2L^-^ + ArtkD,L'a)5qk , (87) 

k=l 



where we have used the fixed end points constraint to derive the last equation 
and defined 



Finally we obtain the modified Euler-Lagrange equations by setting the vari- 
ations to zero: 



ArtkDiLdiqk, - ArD2Ld{qk-i, ^r^^ 



Cfe - ek-i = , 

ArQk-l ^ 

Wo. 
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Lemma 15 The algorithm defined by (88) preserves the Lagrangian two-form 
and the energy. 



PROOF. The first equation of the algorithm proves energy conservation. To 
show that the Lagrangian two-form is preserved, we compute dS^ along a 
discrete trajectory: 



71—1 T-\ jk—l 

dS^ = tY, A^iL'^f'dtk) + A,{D2L',-'dqk) - A,(-^A,gfe_idtfe) 

k=l ^rtk-l 

n—l 



T^Ar{ekdtk + D2L''a ^dqk) 

k=l 

rf^A^e^,, (89) 



k=i 



where 6^ = Ckdtk + D2L^ ^dqk- Since d'^ = 0, we obtain that the symplectic 
two- form Lu^ = d9^ is preserved along the trajectory. 



The proof of this lemma only involves the modified Leibnitz law and does not 

depend on the definition of the discrete Lagrangian function. As a consequence, 
it also applies if one derives modified velocity Verlet and Newmark algorithms. 



5.2.1.2 Hamiltonian approach Let the Lagrangian function be L(g, g) 



\fMq - V(q). Then 



.2Artk A^tk 
and the associated momenta are: 



Pk+i = M——, 

e.+i = --^M^-y(?,). (91) 
2A,tk A,tk ^^"^ ^ ' 

The discrete Hamiltonian function is then: 

Hd = AM^Pl+iM-bk+i + V{qk)) = A^tkHd{qk,Pk+i) ■ (92) 
One can readily verify that Hd{qk,Pk+i) = —^k+i- 
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Let us now derive the modified discrete equations of motion by applying tlie 
GDMHP (tlieorem ( 14) ) . We skip a few steps in the evaluation of the variations 
of to finally find: 



n-l 



6S^ =T6Y,{Pk+l,^Tqk) -Hd{qk,Pk+l) 



k=0 

n-l 



k=0 



-{ArPk + ArtkDiHd{qk:Pk+l):Sqk) + ArCk+lStk+l ■ 

The variations of {6qk,Spk+i,Stk) being independent, we obtain: 



(93) 



A^gfc = ArtkPk+i , 
KPk = -A^tkVV{qk) , 

Arek = 0. (94) 



Lemma 16 The algorithm defined by equations (94) preserves the symplectic 
two-form and the energy. 



PROOF. The proof proceeds as the previous ones, we compute dS^ along a 
discrete trajectory. We skip the detail of the computation: 

n-l 

dS^ ^Tj2^r{Pk,dqk) + ekdtk. (95) 

fe=0 

Define = {pk, dqk) + Ckdtk and = d9j^ . Since = 0, we obtain 

A,u;f = 0. 



Remark 17 The 1-form 6^ corresponds to the contact 1-form 9 encountered 
in continuous time dynamics. Indeed, if one remembers thatck — —Hd{qk-\-iPk), 
then we have: 



9^pdq- Hdt , (96) 
6lf ^Pkdqk - Hd{qk-i,Pk)dtk ■ (97) 
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5.2.2 Midpoint discretization 



In the same manner, we can apply the modified variational principle to other 
discretization. For the midpoint scheme we have — ^*^+^+^'' and the mod- 
ified Leibnitz rule is defined by equation (46). Let us define the generalized 
momenta: 



Pk+i +Pk _ d_ 9Ld 



Pi-^a^ (98) 



Then applying the modified discrete Hamilton's principle (Definition (14)) 
yields (after a few simplifications): 

n-l 

5S^ = rj:{AXD,L', - A^pt, 5qt) - , (100) 

fe=0 



where = ^diQki a^)' variations (5g^, 5i^) being independent, we ob- 
tain: 



Pfe+l — Pk ^fe+l ~ ^fe 7-V T ( + l + ^k + l — Qk V 

T T 2 tk+l - tk 

Pk+1 + Pk tk+l ~ tk J- I Qk+l + Qk Qk+l — Qk x 
9 ~ ^2^dl 5 , —) , 

gfc+i + ^fc _ ^ ^ gfc+1 + gfc gfc+i ~ gfc 
2 2 tfc+i — tk 

-{D,U{^^^^. ?^±^), . (101) 

Lemma 18 The algorithm defined by equations (101) preserves the Lagrangian 
two-form as well as the energy. 



PROOF. We omit the proof since it proceeds as before. 



Now define the discrete Hamiltonian function Hfi{q'l,pf) = iJ(^^^^y^, l^£±i±i^) 
and the modified Hamiltonian function H^i = A'^tfH(i{qk,pi)- Then applying 
the GDMHP yields: 
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Qk+i — Qk tk+i —tk„ „ . Qk+i + Qk Pk+i + Pk X 
r r ^ 2 ' 2 ^ 

Pk+l — Pk tk+l ~tkr^ ^ t + <lk Pk+1 + Pk x 

—r- ^ — ^^i^<^(^^' ' 

Cfc+i — Cfc = , 

2 ^"^"^^ 2 ' 2 

Lemma 19 T/ze algorithm defined by equations (102) preserves the symplectic 
two-form as well as the energy. 



PROOF. We omit the proof since it proceeds as before. 



5.3 Concluding remarks 



The algorithm defined by equations (102) is the same as the one developed 
by Shibberu [32]. Shibberu's approach corresponds to the first formulation of 
the GDMHP (definition (13)) for the midpoint rule but he used a different 
discrete variational principle from DVPII. 

One other work on symplectic energy preserving algorithms is that of Kane, 
Marsden and Ortiz [23] . They developed a generahzed discrete modified Hamil- 
ton's principle that is based on DVPI. Their approach is different from ours: 
they assume a different time step at each iteration, and then take the variation 
of the discrete action without varying the time step (i.e., in a n dimensional 
space). As a consequence they only obtain n equations for the n + 1 variables 
{Qk: hk) where hk is the time step at the /c*'* step. They then add an energy 
constraint to obtain n + 1 equations. Their definition of the energy is similar 
to ours and therefore both methods provide the same algorithms. However, 
there are fundamental differences between the two methods. First, the method 
developed in this paper is fully variational. Second, all the differences between 
DVPI and DVPH that we emphasize at the beginning of this paper still remain 
because their work is based on DVPI whereas our is based on DVPII. 



6 Discrete Hamilton-Jacobi theory 



So far we have developed two variational principles that are the discrete coun- 
terparts of Hamilton's principle on the tangent bundle and on the cotangent 
bundle. Through several examples we have observed that both variational 
principles are equivalent and that they allow us to recover classical variational 
symplectic integrators. We have also shown that they can be modified so that 
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energy conservation is assured. In this section, we concentrate on discrete 
Hamilton- Jacobi theory. We define discrete canonical transformations (DCT), 
discrete generating functions (DGF) and derive a discrete Hamilton- Jacobi 
equation that allows us to show that the energy error for a certain class of 
scheme is invariant under discrete canonical transformations. 



6.1 Discrete symplectic geometry 

We consider again a discretization of the time t into n instants T = {{tk)kG[i,n]} 
but we restrict here to the case where Mk is a n-dimensional vector space. We 
still define M = \JMk. 

Definition 20 A discrete symplectic form uj on Ai is such that at tk, OJ — ujf, 

where is a non degenerate, closed, two-form on Mj} = U Mk-^-i- 

A discrete canonical one-form, on M. is such that at tk, — Of, and uj'l — 

-dot 

A discrete symplectic vector space (A^, cu) is a vector space M. = [j Mk together 
with a discrete symplectic two form on M. . 

Using a symplectic chart, a discrete symplectic form on M. at tk can be written 

as: 



and the canonical one-form as 6*^ = pfdq'jf. 

In the remainder of this section we consider the geometry associated with the 
midpoint scheme, that is, we define z'^ = {qk,Pk) — ^''^^^^^ and use 

the modified Leibnitz law (46). However, the content of this section can be 
applied to any scheme as long as one can define a discrete Hamiltonian vector 
field from the discrete Hamiltonian function and the discrete symplectic two- 
form (see next definition). It is clear that the theory herein can be adapted 
to systems for which the action integral involves a term of the form Hd{zf.), 
where linear combination of Zk and Zk-^-l but it is not clear if it can be 

adapted to the Stormer rule for instance {zf = {qk,Pk+i) cannot be written 
as a linear combination of Zk+i and Zk so the next definition does not apply). 
We do not know how to modify this approach so that a discrete Hamiltonian 
vector field can be defined from the Hamiltonian function Hd{qk,Pk+i)- 

Definition 21 Let (A^,a;) be a discrete symplectic vector space, and Ha : 
Al — > R a smooth function. Define the discrete vector field Xjj such that at 
tk, Xfj — X^, where X^ is of the form 



uji = dqt A dpi , 



(103) 



d 



(104) 



X 



dpi 
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and verifies: 

ix^ui^dHa. (105) 

The discrete vector field Xfj is called the discrete Hamiltonian vector field, 
is called a discrete Hamiltonian system. 

Proposition 22 Using the canonical coordinates, a Hamiltonian vector field 
is of the form: 

Xfj^J-dHa, (106) 



where J = 



\-">/ 



PROOF. Equation (105) is expressed in local coordinates as: 

i^a {dqi A dpi) = D,Ha{zi)dqt + D,H,{zt)dpt . (107) 



Let Xfj be: 



Xj^A?,f^ + AM^, (108) 



then, 



ixf^ {dqt A dpi) = {ixj^dqi)dpl - dqt A {ixj^dpt) (109) 
^AUtdpt-Aiptdqt (110) 

Identifying this last equation with equation (107) leads to equation (106). 
6.2 Discrete canonical transformation 



We now define the class of discrete canonical transformations. The definition 
given here is restricted to linear time-dependent maps (with respect to the 
phase space variables) . We believe larger class of transformations may be con- 
sidered if one works with discretization of the spacetimc [25]. Let (A4i,a;i) 
and (7^2)^2) be discrete symplectic vector spaces and JF be the set maps 
/ : T X Ail — >■ T X that are linear with respect to the phase space vari- 
ables. Consider a map f ^ T such that Vt^ G T, f{tk,-) = fk{-) where fk is 
the following linear map: 



— {Qk,Pk) ^ Zk — {Qk, Pk) — ^k^k + Bk . 
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Since fk is linear, we have: 



fk{zt) = l{fk{zk) + fk{zk+i)), (111) 



/,(A^4) = AA^4. (112) 

Definition 23 A linear, time- dependent map f is called a discrete canonical 
transformation (DCT) (or a discrete symplectic map) if and only if f*uj2 — oJi, 
or equivalently, V/c e [l,n], fk^^2,k — ^i,k- 

Proposition 24 If f is a DCT then is invertible for all /c e [1, n] 



PROOF. Suppose there exists a k such that Ak is not invertible. Then 3z^ e 
such that 

3v,eT,,Ml,\Tfk-v,=0. 

Then, \/v2 G T^dML\v2 7^ 0, ujf Jvi,V2) = u^iki^fk ■ vi,Tfk ■ V2) since / is 
symplectic. The right hand side is zero but the left hand side is not. This is a 
contradiction and therefore Ak is invertible. 

Lemma 25 Let f be a discrete canonical transformation. Then fk^^2k — ^ik 
can be written in the matrix form AkJA]^ = J. In addition, f preserves the 
form of the discrete Hamilton's equations. 



PROOF. The statement AkJA]^ = J is just the matrix statement of /fcCjf ^ = 
cuf j^.. Let us prove that / preserves the form of the discrete Hamilton's equa- 
tions. Define the function such that K^o f = H^. 

On one hand, using equation (112) we have: 



\d _ fk{Zk+l) fkjZk) (113) 

^Ak^Ut- (114) 
On the other hand: 



JVHa{zi) = JV{K,of^{zi)) (115) 

^JAlVKaizi). (116) 

Since AkJA^ — J, we obtain: 

Afz', = JVK,{zt) (117) 
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This last result can be summarized as follows: 

Proposition 26 Let Xfj be a discrete Hamiltonian vector field with Hamil- 
tonian function and f a discrete symplectic map. Then f^Xfj is a discrete 
Hamiltonian vector field with Hamiltonian function f*Hd. 



6.3 Discrete generating functions 



Proposition 27 Let (A^i,6<;i) and (Al2,<^2) &e two discrete symplectic vector 
spaces, TTj : M.i x M.2 M-i the projection onto Aii and define 

VL = TtIuJi — TT2W2 ■ (118) 

Then, 

(1) Q is a discrete symplectic form on Aii x Ai2, 

(2) a map / : A^i — >• M.2 is a discrete symplectic map if and only ifi*jrfl = 0, 
where if : T f —>■ Mi x M.2 is the inclusion map and Tf is the graph of f. 



PROOF. We recall that at t^, Q = where = ^lk^fk-^2k^2k- To prove 
that Q is a discrete symplectic form, we need to prove that is a symplectic 
form on Mf ,^ x M^^. for all k e[l,n]. 



dnt^dirrlujl^-Tr^ul,) (119) 
^nldu;l,-n*2dui, (120) 
= 0, (121) 

since cuf^j^ is closed and d commutes with the pull back operator. 

Now let zt = «fc,4,) e Ml, X and v = {vi,V2) G T^^iM^k x M^,) - 
T. Mf, X T^d M^, such that 

yw = {wi, W2) e 7;.(M^, X Ml,) nt{v, w) = o (122) 

and let us prove that v is zero. We have 



nt{v, w)^u;l,{m{zt)){Tm ■ V, Tm ■ w) - u;l,{7r2{zt)){Tn2 ■ v, Tn2 {m 
-<ki4,k)ivi,wi)-u;l,izl,)iv2,W2) (124) 
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The right hand side of equation (124) is zero for all w if and only if both terms 
are zero, that is, 

u;l,{zl,){v^, wr)^0, <fc(4fc)(^^2, w^) = (125) 
Since cuf^^. is non degenerate, vi = V2 = and flf is closed. 

We now prove the second statement of the proposition. We first notice that 
fk induces a diffeomorphism of M^ p. to F/^, , so we can write 

Ti4,M4}) = {(^'^/^ • e T^i^l,} (126) 

Then, 

entiivi, Tfk ■ V,), {V2, Tfk ■ V2)) =<fc(^l, V2) - uikiTfk ■ Vl, Tfk ■ V2) 

-{<k-fKk){vuV2) (127) 

Hence, fk is symplectic if and only if = 0, i.e., / is a discrete symplectic 
map if and only if i*Q, — 0. 

Using the Poincare lemma we may write = —dQ'l and the previous propo- 
sition says that i^^t is closed if and only if / is a discrete symplectic map. 
Using again the Poincare lemma, we conchide that if / is a discrete symplec- 
tic map then there exists a function S : T f ^ M. such that i^O = dS, i.e.. 

Definition 28 Such a function S is called a discrete generating function for 
the discrete symplectic map f.Sis locally defined and depends on the choice 
ofO. 

• Let elk = Pkdqt and = P^dQi then 

dS = —{ql Qt)dqt + ^{qt Qt)dQt , (129) 

that is, 

(99 (99 
Pt = Q^iQt, Qt) Pt = Qi) ■ (130) 

S as defined corresponds to a discrete generating function of the first kind. 

• Let 01, = pfdqi and 0^, = -QfdP^, then 

t}^ei = (TTi o tfjY.dqt + (vr^ o Zf.YQidP^ , (131) 

(99 (99 
dS-g^iqt, Qt)dqt + g^iqt, Q'M , (132) 
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that is, 

pt = ^i<iiQt) Qt-%{<iiPt)- (133) 

S as defined corresponds to a discrete generating function of tlie second 
kind. 

In tlie same way, one can define 4" generating functions as in tlie continuous 
case. Note tliat since / is linear witli respect to its spatial variables, S is 
also linear with respect to its spatial variables. At tk, S = Sk where Sk{-) ~ 
Tki') + Uk is affine map, is a 2n x 2n matrix and C/fe is a 2n x 1 matrix. 



6.4 Discrete Hamilton- J acohi theory 

In this section we use the notions introduced previously to develop a discrete 
Hamilton- Jacobi theory. Let / be a discrete symplectic map, let Mff^ = T*Qfi^ 
and let S be an associated discrete generating function such that ai tk, S = 
where Sk{-) = Tk{-) + Uk 

Theorem 29 Define 

Ptiqi Qt) = D,Sk{qt Qi) , PM, Qi) = -D2Sk{qt Qt) ■ 
Then the following two conditions are equivalent: 

(1) S is a discrete generating function associated with f; 

(2) • For every curve {ck)k in Qi — U Qi,fe satisfying: 

Aici = Tn*Q,Xf,{ci,pt), (134) 

the curve k ^ (Cfc,^^) is a discrete integral curve of Xfj, where yr^d is 
the cotangent bundle projection onto the configuration space. 
• For every curve {ck)k in Q2 — {j Q2,k satisfying: 

Ar4 = Tn*Q.Xi{ciPt), (135) 

the curve k 1— > {cf,P^) is a discrete integral curve of Xf., where tTq^ 
is the cotangent bundle projection onto the configuration space. 



PROOF. Suppose 5 is a discrete generating function, let Q'^ be fixed and 
consider a curve {ck)k verifying 

Aict^T7r*Q,^XUct,Pt), (136) 
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In other words, verifies: 

Aici = D,H{4,pt), (137) 

Since 5" is a generating function, pf is the momentum associated with cf and 
verifies: 

Aipt = -D,H{ct,pt). (138) 

These last two equations are exactly a restatement of: k i-^ {ctiPk) ^ dis- 
crete integral curve of Xfj. To derive the second item we proceed in the same 
manner, but this time is fixed. 

Now we suppose item (2) and we show that 5 is a discrete generating function 
for /. The statements k h- > (cf,p^) is a discrete integral curve of Xfj and 
A; I— > (c^, P^) is a discrete integral curve of Xf^ are equivalent to saying that 
p'l and are the momenta associated with the generahzed coordinates, and 
therefore, 5" is a generating function for /. 

Theorem 30 We consider again a time dependent function S which is linear 
with respect to the spatial variables. Then the following two statements are 
equivalent: 

(1) S is a discrete generating function associated with f; 

(2) For every H there is a function K such that 

H{qi D,S{ql QD) = K{Qi D,S{ql Qf)) (139) 

PROOF. Suppose is a discrete generating function. Then from the previous 
theorem, for every curve (c^, Ck) in Qi x Q2 satisfying Atci = Tt!-*^ Xfr(ci,pi) 

and AfC^ = Ttt* , X|.(C^, P^), the curves k ^ (cipi) and k ^ (Cf , P^) are 

discrete integral curves of Xfj and respectively. Then, using the symplectic 
identity ([1] page 382) that holds for any function S 

ul^D^S o n*Q,J -v^w)^ ui^{v, w - T{D,S o n*Q,J ■ w) 

we get: 

ul,{T ( D,S o n*Q,J . Xfj{ck, DiSk),w) = 

<fe ( Xfj(ck, D^Sk),w) - dHaick, D,Sk) ■ TD,S{ck, D,Sk) ■ w (140) 
u;l,{T ( -D2S on*Qa)- X|-(Cfc, -D2Sk),w) = 

(141) 
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In addition, since pi = -Di5(cf , C^) and = -DiS{ci, C(), 



Aipt = TD,S{ci Ct)Aici = T{D,S o tt^, ) • X^(c,, D,Sk) (142) 
A^P,'' = T(-Z}2^o7r^, )-Xi(Cfc,-D2^fc). (143) 

f being a discrete canonical map, Tf^i/^fpf) — ^fP^ so the left hand side 
of equation 141 is the image under f of the left hand side of (141). Using 
proposition (26), we conclude that: 

Tfk-dHd(ck, D,Sk)-TD,S(ck, D^Sk) = -dKa(Ck, -D2Sk)-TD2S(Ck, -D^Sk) , 
which is equivalent to the discrete Hamilton- Jacobi equation. 
The proof that 2. implies 1. follows from these arguments. 

6.5 Applications of the discrete Hamilton- Jacobi theory 

The goal of this section is to highlight the benefit of having a discrete Hamilton- 
Jacobi theory. First, we have proven the invariance of the discrete Hamilton's 
equations under a certain class of coordinate transformations. Second, we have 
shown in theorem 30 that changing coordinates using a discrete symplectic 
map does not improve the performance of the algorithm in terms of energy 
conservation. As a consequence we have the following lemma: 

Lemma 31 The midpoint scheme preserves the energy for linear systems. 

PROOF. The discrete phase flow for linear systems is piecewise hnear contin- 
uous and the map {qk,Pk) ^ iQk+i,Pk+i) is symplectic (the midpoint scheme 
is a symplectic algorithm). Therefore, the discrete phase flow is a discrete sym- 
plectic map that maps H into a constant K. Integration of the new Hamilto- 
nian system defined by K is trivial {{Qk+i, Pk+i) = {Qk,Pk)) and obviously 
preserves the energy. As a consequence the integration of the Hamiltonian 
system defined by H also preserves the energy. 

Finally, we illustrate the use of the above material with a nonlinear example. 
We study the energy error in the integration of the equations of motion of a 
particle in a double well potential using different sets of canonical coordinates. 
Consider a particle in a double well potential, i.e., H = + |(9^ ~ 
As shown in figure (2), the midpoint scheme does not preserve the energy. 
The following time-dependent discrete canonical transformation (at each step 
the transformation is a different expression) Zk — AkZk + Bk where Ak — 
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/ cos(A;^) -sin(fc^)\ 

. and Bk — rotates the system by k9 — /carccosO.99 at 

\^sin(A;^) cos(A;^) J 

the /c*^ step. In figure (3) we plot the same trajectory in the new system of 
coordinates, the energy error is exactly the same. In other words, the energy 
error is invariant under discrete canonical maps. 
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(a) Trajectory of in the q—p plane (b) Energy error for constant 

time step midpoint scheme as a 
function of time. 



Fig. 2. Particle in a double well potential with initial conditions {q,p) = (1,0.05) 
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(b) Energy error for constant 
time step midpoint scheme as a 
function of time. 



Fig. 3. Particle in the Hamiltonian vector field f*Xfj, where Xfj is the Hamilto- 
nian vector field corresponding to a double well potential. Initial conditions are 
(Q,^')=/o(l,0.05). 



7 Optimal control 



For a general optimal control problem, necessary conditions for optimality 
may be derived from the Pontryagin maximum principle. These conditions 
often yield equations of the same form as Hamilton's equations coupled with 
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nonlinear equations. We have seen previously that Hamiltonian systems, i.e., 
Hamilton's equations, can be integrated using symplectic integrators. How- 
ever, if Hamilton's equations are coupled with algebraic nonlinear equations, 
the above theory does not apply. What is the correct discretization of the alge- 
braic equation? In this section, we develop a discrete maximum principle that 
tackles this problem and provides a unified view on solving optimal control 
problems using symplectic integrators. 



7. 1 Necessary conditions for optimality 



7.1.1 Problem Statement 

Let J = g{x,u)dt be a performance index (also called a cost function) and 
consider the following optimal control problem: 



min J g{x,u)dt, (144) 
subject to the dynamics 

x^f{x,u), (145) 
and to the initial and final time constraints: 

(l>Mto):to) = 0, 0/(a;(t/),t/) = 0, (146) 

where / and g are functions from x to M of class C^. 



7.2 Maximum principle 



To solve the optimal control problem, we apply the maximum principle. 

Theorem 32 (Mciximum principle) Solutions to the optimal control prob- 
lem defined by equations (144), (145) and (146) correspond to critical points 
of the cost function J in the class of curves 7 = {x{t),u{t)) G F where F is 
the set of curves satisfying (145) and (I46). 

Remcirk 33 This formulation differs from the one given by Pontryagin [30] 
but the main point of the Pontryagin maximum principle is that it yields neces- 
sary conditions for optimality under far less severe regularity conditions. The 
above formulation is based on the equivalence between the Pontryagin max- 
imum principle and the calculus of variations in the case where the control 
region is an open set in a finite dimensional vector space (see [30] chapter V 
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for more details). It is therefore equivalent to classical variational formulations 
given in Block et al. [3,4] and Gregory and Lin [11] for instance. 



To apply the maximum principle we first need to define the augmented cost 



function J^: 



Ja= 9{x, u) + {p,x- f{x, u))dt + (Ai, <pi{x{to),to)) + (A/, (j)f{x{tf),tf)) 

J to 

= / H{x,p,u) - {p, x)dt + {Xi, (f)i{x{to),to)) + {Xf,(i)f{x{tf),tf)) , 

J 

where the p's, the Aj's and the A/'s arc Lagrange multipliers and H[x,p,u) — 
g{x, u) + (p, /(x, u)). Taking variations of the augmented cost function assum- 
ing fixed initial and final time yields: 

SJa^S (^j^ ' H{x,p,u) + {p,x)dt^ + S{Xi,(f)i{x{tQ),to)) 
+5{Xf,ct>f{x{tf),tf)) 
= / {D2H{x, p, u) - X, Sp) + {DiH{x, p, u) + p, 5x) 

Jto 

+ {D^H{x,p,u),5u)dt+ {-p{tf) + Di4PjXf,5xf) 

+ {p{ti) + D^(t)JXi, 5xi) . 

We now let the variations of Ja be zero to obtain necessary conditions for 
optimality: 

x = D2H{x,p,u) , (147) 
p=-DiH{x,p,u) , (148) 
0^D3H{x,p,u), (149) 

as well as transversality conditions: 

p{ti) = -Di<Pi{x{to),tofXi , p{tf) = Di<Pf{x{tf),tffXf . (150) 

Equations (147)- (150) define the necessary conditions for optimality. 

7.3 Solving the necessary conditions for optimality 



To solve these conditions, the most common technique is to find the optimal 
control feedback law from (149) and then use a shooting method to solve the 
two-point boundary value problem defined by (147), (148) and (150). More 
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precisely, suppose (149) allows one to solve for u as a function of {x,p) and 
define the Hamiltonian function 



H{x,p) ^ H{x,p,u{x,p)) , (151) 
then the necessary conditions (147) and (148) simplify to: 

x = D2H{x,p), (152) 

p=-DiH{x,p). (153) 



Equations (152) and (153) define a Hamiltonian system that has no physical 
meaning in general. As we will see later, for sub-Riemannian optimal control 
problems the Legendre transform is ill-defined and therefore DVPI cannot be 
used to discretize such systems whereas one could use DVPII (theorem 3). 
However, one may not be able to solve (149), and then the question of how 
one can use symplcctic integrators to solve the optimal control problem arises. 
What is the correct discretization of (149)? In the next section wc address this 
issue. Specifically, we introduce a discrete maximum principle that allows us 
to derive discrete necessary conditions for optimality that are in agreement 
with the one obtained from the maximum principle. 

7.4 Discrete maximum principle 

7.4-1 Problem statement 

In discrete settings, the cost function is 

n-1 
k=0 

and the optimal control problem (144) is formulated as: 

n-1 

minY^ gd{xi,ui)T, (154) 

"fc k=0 

subject to the dynamics 

Afxt^Uxiut), (155) 

and to boundary conditions: 

where fd and ga are functions from R" x R"^ to R of class C^. They correspond 
to discretization of the continuous time functions / and g. 
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H.2 Discrete maximum principle 

To obtain necessary conditions for optimality, we define the following discrete 
maximum principle, the discrete counterpart of the maximum principle: 

Definition 34 (Discrete mclximum principle) Solutions to the discrete op- 
timal control problem correspond to critical points of the cost function J in 
the class of discrete curves 7 G F, where F is the set of all discrete curves 
ixk,Uk)keli,n] that verify (155) and (156). 

Remark 35 The above definition is the discrete counterpart of the maximum 
principle. It compares to previous works on discrete optimal control theory that 
extend the Pontryagin maximum principle to discrete systems such as Jordan 
and Polak [22] as theorem 32 compares to the Pontryagin maximum principle. 
In other words, in contrast with Jordan and Polak [22], we restrict the class of 
discrete optimal control problems so that we can derive necessary conditions 
that define symplectic algorithms. 

As in the continuous case, to find critical points of J under the non-holonomic 
constraint defined by equation (155), we must append the constraints to J 
using the Lagrange multipliers. The resulting function is called the augmented 
cost function: 



Ja = T.(9d{xt, 4) - {pi Afxt - f.ixi, ut)))r + (Ao, 0o) + (A„, (157) 

fc=0 

j£{Ha{xipiui) - (p^, A^x^))t + (Ao,0o) + {KAn) , (158) 

k=0 

where the PkS, the Aq's and the A„'s are Lagrange multipliers and Hfi{xf, pf, uf) = 
gd{xf,uf) + (pf., fd{xf,uf.)) . To apply the discrete maximum principle, one 
needs to specify the discrete derivative operator as well as the expressions of 
xf., uf and pf as a function of {xk+i, Xk), {uk+i, Uk) and {pk+i,Pk) respectively. 



7.4.3 Examples 

7.4.3.1 Stormer's rule If we choose to be the forward difference A^- 
and {xf,pf,uf) = {xk,Pk+i,Uk) then we recover the discrete maximum princi- 
ple developed by Bloch, Crouch, Marsden and Ratiu [5]. 
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\fe=0 / 

n-1 

= ^{D2Hd{Xk-,Pk+l-,Uk) - ArXk,Spk+l)T 
k=0 

+ {DiHd{xk,Pk+i,Uk) + A^pk,Sxk)T{+DsHd{xk,Pk+i,Uk),Suk)T 

+ (00, (^Ao) + {(f>n, SXn) + {-Pn + -Di0^A„, SXn) + {po + -Di0o Aq, Sxq) , 

(159) 

where the modified Leibnitz law (1) has been used. We impose the variation of 
the augmented cost function to be zero to obtain discrete necessary conditions 
for optimahty and transversahty conditions: 



A^Xk = D2Hd{xk, Pk+i, Uk) , (160) 

ArPk = -DiHd{xk,Pk+i,Uk) , (161) 

= D3Hd{xk,Pk+i,Uk) , (162) 

Po = -Di(j)o{xo, to)'^Xo , Pn = Di(j)n{Xn, t^f K ■ (163) 



The algorithm defined by (160), (161) and (162) is equivalent to the one 
derived by Bloch, Crouch, Marsden and Ratiu [5] for the symmetric rigid 
body. 

Lemma 36 The algorithm defined by (160), (161) and (162) is symplectic. 
PROOF. Define the cost function as: 

n— 1 

Ja = J2i^d{Xk,Pk+l,Uk) + ipk+U A^Xk))T . (164) 
fe=0 

Ja is the augmented cost function from which we have removed the boundary 
conditions. Boundary conditions yield transversahty conditions, that is con- 
ditions on the initial and final states of the system. Hence these terms are 
irrelevant to the study of the advance map {xk,Pk,Uk) ^ {xk+i,Pk+i,Uk+i). 
As in discrete dynamics, we consider d'^Ja, assuming {xk,Pk,Uk) verifies the 
above necessary conditions and we obtain: 

n-1 

dJa = ^ A-r(pfe, dxk)T . (165) 

fe=0 

Prom = 0, we conclude: 

n-1 

= ^ Ardipk, dxk)T , that is, V/c e [0, n — 1] , dpk-^-i A dxk+i — dpk A dxk ■ 

k=0 

(166) 
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The symplectic nature of the algorithm is obtained directly from the varia- 
tional principle - there is no need to compute dpk A dx^ and dpk+i A dx^+i- 



7.4.3.2 Midpoint scheme Midpoint discretization may also be obtained 
if we choose 

^ Xa + i + Xk a Pk+1 +Pk d '-''A+l + Uk 

•^fc ~ 2 ^ J^k — 2 ' ~ 2 

and = Rr/2 — R-t/2- One can readily verify that the discrete maximum 
principle yields the following necessary conditions for optimality and transver- 
sality conditions: 



Aixt = DMxt,pt,ut), 
Aipi = -D,H4xtpi,4), 
O^D,Hd{xtptut), 

PO = Di(j)o{xo, to)^Xo , Pn = --Dl0„(x„, tnfXn ■ 



(167) 
(168) 
(169) 
(170) 



Lemma 37 The algorithm defined by (167), (168) and (169) is symplectic. 



PROOF. We omit the proof since it proceeds as before. 



7.5 Discrete maximum principle v.s. discretization of the Pontryagin maxi- 
mum principle 



So far we have considered two methods for obtaining a symplectic algorithm 
that integrates the necessary conditions for optimality. The first method, which 
applies only to a certain class of problems, consists of discretizing the necessary 
conditions obtained from the Pontryagin maximum principle once the control 
as been expressed as function of {x,p). The second method consists in using 
the discrete maximum principle. In this section, we show that under certain 
assumptions both methods are equivalent, that is we prove the commutative 
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diagram (171). 



X = f{x,u) 



{PMP) 



{BMP) 



(171) 



H{x,p,u) (DMHP) Hd{xi,pi,uf) 



■> 



H{x,p) Hd{xi,pi) 



where H is defined by (151), DMHP stands for discrete modified Hamilton's 
principle, PMP stands for Pontryagin maximum principle, and DMP stands 
for discrete maximum principle. 

We recall the required assumptions to prove the equivalence of the diagram. 
We assume that (149) can be solved for as a function of {x,p) and that the 
initial and final states x{tf) — Xf and x{to) — Xi are given. In addition, we 
impose ga ^ g and fa ^ f . 

To discretize the Hamiltonian system defined by H, we use the discrete mod- 
ified Hamilton's principle: 



for any variations of {xf.jpf) and Sxq = Sxn = 0. One can readily check that 
(172) can also be written in an equivalent form as: 



for any variations of {xf.,pf., uf) and Sxq = Sxn = where uf is now considered 
as an independent variable. In addition since f = fd and g = ga, H = Ha, and 
we conclude that the discrete modified Hamilton's principle as formulated and 
the discrete maximum principle are equivalent. 

7. 6 The Heisenberg optimal control problem 

The Heisenberg problem (Brockett [6], Bloch et al. [3]) refers to under actu- 
ated optimal control problems which are controllable. For instance, consider 
a particle that has two actuators in the {x, y)-p\ane and with velocity in the 




(172) 




(173) 
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z direction defined by i = — xy. This system is controllable, however, to 
reach a point (a > 0, 0, 0) from the origin (0, 0, 0) requires a non-trivial con- 
trol vector. In the following, we study the Heisenberg problem to illustrate the 
approaches we have developed above. This problem formulates as: 



min f\u,u)dt, (174) 

U=(ui,U2) J to 



subject to 



x — u, (175) 
y^v, (176) 
z — uy — vx, (177) 

and to the boundary conditions: 

{x{to),y{to), z{to)) = (0, 0, 0) , {x{tf),y{tf), z{tf)) = (a > 0, 0, 0) . 

This is a hard constraint problem, therefore the transversality conditions are 
of no use; They yield 2n equations but introduce 2n new variables. 

Define if as ^ 

H{q,p,u) = -{u,u) + {p,q) , 

where q = {x,y,z) and p = {Px,Py,Pz)- The Pontryagin maximum principle 
yields: 



dH , , , , 

q^^{q,p,u), (178) 

OH 

P = —g^{Q,P,u) , (179) 
dH 

0^—{q,p,u). (180) 
Equation (180) allows us to solve for u as a, function of (g,p): 



Ul^Px+ PzV , U2^Py- PzX , (181) 

Hence, equations (178)- (179) become: 



^ = (182) 
dH 

P^—g^iq,p), (183) 
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where 



H{q,p)^H{q,p,u{q,p)) 
1 



- O (Px + Py) - PxPzV + PyPzX 



(184) 



Equations (182) and (183) are of the same form as the Hamilton equations. 

Therefore, the necessary conditions for optimahty yield a Hamiltonian system 
with Hamiltonian function H. We now prove that H is degenerate at the 
origin, and so is the Legendre transform. The Hessian of H is: 



' dH ' 



/ 



-1 
-1 X 
-y X 



y 



Thus, det 



dH 

d{q,p) 



X^ + y^. 



i.e., the determinant of the Hessian of H is 



singular at (0, 0). As a result, it is not, a priori^ possible to define a Lagrangian 
function associated with the Hamiltonian H using the Legendre transform ^ . 
Therefore, the discrete modified Hamilton's principles (DMHP) must be used 
to discretize Eqns. (182) and (183). One cannot use a discrete Hamilton's 
principles (DHP) for instance because the system is not Lagrangian. This point 
is of importance. It motivates the need to introduce the variational principles 
presented in this paper, as previous works on variational principles mostly 
focused on systems with non-degenerate Lagrangian functions. To discretize 
the necessary conditions, we choose the geometry associated with the Stormer 
rule and using DMHP (definition 3) to eventually find the following symplectic 
algorithm: 



Arqk^D2H{qk,Pk+i) , 
ArPk^-DiH{qk,Pk+i) ■ 



(185) 
(186) 



Let us now discretize the Heisenberg problem using the second approach, based 
on the use of the discrete maximum principle. We first discretize the problem 
statement: 

-j^ n— 1 
Wfc=("l,fc,'"2,fc) ^ 

subject to 



(187) 



° Using Lagrange multipliers one can define a Legendre transform and find a La- 
grangian function associated with the system. We refer to Bloch [3] for a presentation 
of this technique that involves variational principles with constraints. 
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ArXk 
^rVk 
^rZk 



= U2,k , 
= Ul,kyk 



- U2,kXk ■ 



(188) 
(189) 
(190) 



Define the discrete augmented cost function J^: 

n-l 

Ja = Y,^d{qk,Pk+l,Uk) - {pk+U^rqk) , (191) 
fe=0 

where Hd{qk:Pk+i:Uk) = {uk,Uk) + {Pk+uQk) and = {xk.Vk.Zk). To find 
discrete necessary conditions for optimahty we set the variations of J a to zero, 
and we obtain: 



Ar^fe = D2Hd{qk, Pk+i, Uk) , (192) 

ArPk = -DiHd{qk,Pk+i,Uk) , (193) 

Q^D3Hd{qk,Pk+i,Uk). (194) 

Equation (192) allows us to find Uk as a function of {qk,Pk+i)' 

Ul,k = Px,k+1 + Pz,k+\yk , U2,k = Py,k+1 - Pz,k+\Xk ■ (195) 

We then substitute these expressions into equations (192)-(193): 

Arqk = D2Hdiqk,Pk+i) , (196) 

^rPk = -DiHd{qk:Pk+i) , (197) 



where Hd{qk,Pk+i) = Hd{qk,Pk+i,Uk{qk,Pk+i))- By virtue of the commutative 
diagram, (196) and (197) define the same symplectic algorithm as (182) and 
(183). 

In this example, we chose a trivial discretization of the dynamics and of the 
cost function; f = fd and g = Qd- Other algorithms may be obtained using 
nontrivial discretizations. In that case the equivalence principle may not hold 
but the algorithm we obtain will still be symplectic. In addition, in this ex- 
ample we did not take into account any boundary conditions since we have 
seen earlier in the paper that both methods yield comparable transversality 
conditions. Finally, as in discrete dynamics, the discrete maximum principle 
may be modified in order to yield symplectic-energy conserving algorithms. 
We add an independent parameter and consider the time as a generalized 
coordinate, the optimal control problem then formulates as follows: 

n— 1 n— 1 

miriu Yl 9d{xt ui) {tk+i - tk) = 9d{xt ui)ArtkT . (198) 

k=0 k=0 
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subject to the dynamics 

Aixt = AittU{xt,ut) 



(199) 



8 Conclusions 

In this paper we have presented a general framework that allows one to study 
discrete systems. We have introduced variational principles on the tangent and 
cotangent bundles that are the discrete counterpart of the known principles of 
critical action for Lagrangian and Hamiltonian dynamical systems. We have 
shown that they allowed us to recover most of the classical symplectic algo- 
rithms. In the future, we will try to derive additional symplectic algorithms 
such as the symplectic partitioned Runge-Kutta algorithm. In addition, we 
have seen that by increasing the dimensionality of the configuration space, 
symplectic algorithms may be transformed into symplectic-energy conserv- 
ing algorithms. When time is a generalized coordinate, the dynamical system 
is subject to an energy constraint and we are able to adapt our variational 
principles to take into account such a constraint. In the same manner, our 
approach may be modified to derive symplectic algorithms to integrate non- 
autonomous dynamical and control systems with (non-holonomic) constraints. 
We have also identified a class of coordinate transformations that leaves the 
variational principles presented in this paper invariant and developed a dis- 
crete Hamilton- Jacobi theory. This theory allows us to relate the energy er- 
ror in the integration using different set of coordinates. Finally, for optimal 
control problems we have developed a discrete maximum principle that yields 
discrete necessary conditions for optimality. These conditions are in agreement 
with the usual conditions obtained from Pontryagin maximum principle. In 
future research, we want to use the general framework introduced in this pa- 
per to develop variational principles for multi-symplectic algorithms, that is a 
spacetime discretization will be used instead of the time discretization. Such 
a formulation would allows us to develop efficient numerical algorithms for 
simulation of the motion of rigid bodies and complex interconnected systems. 

Acknowledgement: We would like to thank Jerry Marsden for valuable dis- 
cussions. 
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